# https://chat.deepseek.com/a/chat/s/ea1ec431-f27c-4b80-866c-b355fd7596a6
# https://chat.deepseek.com/a/chat/s/f22decc5-062b-41a3-92f1-44e6e33d0c63
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split, cross_val_score, StratifiedKFold
from sklearn.linear_model import LogisticRegression
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import (accuracy_score, confusion_matrix, classification_report,
                             roc_auc_score, roc_curve, precision_recall_curve, average_precision_score,
                             precision_score, recall_score, f1_score)
from sklearn.feature_selection import SelectKBest, f_classif
import warnings

# 忽略警告
warnings.filterwarnings("ignore")

plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']  # 指定默认字体
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题


# 设置随机种子以保证结果可复现
np.random.seed(42)

# 1. 生成模拟数据集
print("步骤1: 生成模拟数据集")
n_samples = 368
features = [
    'rectal_temp', 'pulse', 'respiratory_rate', 'pain_level', 'abdominal_distension',
    'nasogastric_reflux', 'rectal_mucosal_color', 'capillary_refill', 'peristalsis',
    'abdominal_sounds', 'surgical_lesion', 'lesion_type', 'lesion_location', 'lesion_severity',
    'packed_cell_volume', 'total_protein', 'mucous_membranes', 'surgical_exploration',
    'postoperative_compli', 'antibiotic_treatment', 'analgesics_admin', 'age', 'hospital_days',
    'preoperative_med', 'abdominal_palpation', 'peripheral_pulse', 'endotoxin_level', 'surgery_duration'
]

# 创建空DataFrame
df = pd.DataFrame(columns=features + ['outcome'])

# 生成数据
for i in range(n_samples):
    # 基本生理指标
    df.loc[i, 'rectal_temp'] = np.random.normal(38.0, 0.5)
    df.loc[i, 'pulse'] = np.random.randint(40, 80) if np.random.random() > 0.05 else np.nan
    df.loc[i, 'respiratory_rate'] = np.random.randint(10, 25)
    
    # 疼痛和症状指标
    df.loc[i, 'pain_level'] = np.random.randint(0, 11)
    df.loc[i, 'abdominal_distension'] = np.random.randint(0, 4) if np.random.random() > 0.07 else np.nan
    df.loc[i, 'nasogastric_reflux'] = round(np.random.uniform(0.0, 3.0), 1)
    
    # 临床检查指标
    df.loc[i, 'rectal_mucosal_color'] = np.random.randint(1, 5)
    df.loc[i, 'capillary_refill'] = round(np.random.uniform(1.0, 5.0), 1)
    df.loc[i, 'peristalsis'] = np.random.randint(0, 4)
    df.loc[i, 'abdominal_sounds'] = np.random.randint(0, 3) if np.random.random() > 0.05 else np.nan
    
    # 病变相关指标
    df.loc[i, 'surgical_lesion'] = np.random.randint(0, 2)
    df.loc[i, 'lesion_type'] = np.random.randint(1, 6)
    df.loc[i, 'lesion_location'] = np.random.randint(1, 4)
    df.loc[i, 'lesion_severity'] = np.random.randint(1, 5)
    
    # 实验室指标
    df.loc[i, 'packed_cell_volume'] = np.random.randint(25, 55)
    df.loc[i, 'total_protein'] = round(np.random.uniform(4.0, 8.5), 1)
    df.loc[i, 'mucous_membranes'] = np.random.randint(1, 5)
    
    # 治疗相关指标
    df.loc[i, 'surgical_exploration'] = np.random.randint(0, 2)
    df.loc[i, 'postoperative_compli'] = np.random.randint(0, 4)
    df.loc[i, 'antibiotic_treatment'] = np.random.randint(0, 2)
    df.loc[i, 'analgesics_admin'] = np.random.randint(0, 2)
    
    # 其他指标
    df.loc[i, 'age'] = np.random.randint(1, 20)
    df.loc[i, 'hospital_days'] = np.random.randint(1, 15)
    df.loc[i, 'preoperative_med'] = np.random.randint(0, 3)
    df.loc[i, 'abdominal_palpation'] = np.random.randint(0, 4)
    df.loc[i, 'peripheral_pulse'] = np.random.randint(0, 3)
    df.loc[i, 'endotoxin_level'] = np.random.randint(0, 5)
    df.loc[i, 'surgery_duration'] = np.random.randint(30, 180) if np.random.random() > 0.06 else np.nan
    
    # 生成目标变量
    survival_prob = 0.5
    if df.loc[i, 'pain_level'] > 7: survival_prob -= 0.3
    if df.loc[i, 'packed_cell_volume'] < 30: survival_prob -= 0.2
    if df.loc[i, 'age'] < 5 or df.loc[i, 'age'] > 15: survival_prob -= 0.1
    if df.loc[i, 'surgical_lesion'] == 1: survival_prob += 0.15
    if df.loc[i, 'antibiotic_treatment'] == 1: survival_prob += 0.1
    survival_prob = max(0.1, min(0.9, survival_prob))
    df.loc[i, 'outcome'] = 1 if np.random.random() < survival_prob else 0

# 添加额外缺失值
for col in features:
    df.loc[np.random.choice(df.index, size=int(len(df)*0.05), replace=False), col] = np.nan


# 确保目标变量是整数类型（修复错误的关键）
df['outcome'] = df['outcome'].astype(int)


print(f"数据集维度: {df.shape} (样本数×特征数)")
print(f"目标变量分布:\n死亡: {(df['outcome'] == 0).sum()} ({100*(df['outcome'] == 0).mean():.1f}%)")
print(f"存活: {(df['outcome'] == 1).sum()} ({100*(df['outcome'] == 1).mean():.1f}%)")

# 2. 数据预处理
print("\n步骤2: 数据预处理")
# 分离特征和目标变量
X = df.iloc[:, :-1]  # 所有特征
y = df['outcome']   # 目标变量

# 处理缺失值 - 使用中位数填充
print("  处理缺失值...")
imputer = SimpleImputer(strategy='median')
X_imputed = pd.DataFrame(imputer.fit_transform(X), columns=X.columns)
print(X_imputed)

# 标准化特征 - 逻辑回归对特征尺度敏感
print("  标准化特征...")
scaler = StandardScaler()
X_scaled = pd.DataFrame(scaler.fit_transform(X_imputed), columns=X.columns)
print(X_scaled)

# 3. 特征选择
print("\n步骤3: 特征选择")
# 使用ANOVA F-value选择最重要的10个特征
selector = SelectKBest(f_classif, k=10)
X_selected = selector.fit_transform(X_scaled, y)

# 获取选中的特征名称
selected_features = X.columns[selector.get_support()]
print("所有特征的F值:", selector.scores_)
print("P值:", selector.pvalues_)
print(f"  选中的特征: {list(selected_features)}")

# 4. 划分训练集和测试集
print("\n步骤4: 划分训练集和测试集")
X_train, X_test, y_train, y_test = train_test_split(
    X_selected, y, test_size=0.2, random_state=42, stratify=y
)
print(f"  训练集大小: {X_train.shape[0]} 样本")
print(f"  测试集大小: {X_test.shape[0]} 样本")

# 5. 创建逻辑回归模型
print("\n步骤5: 创建逻辑回归模型")
model = LogisticRegression(
    penalty='l2',           # L2正则化防止过拟合
    C=0.1,                  # 正则化强度
    solver='lbfgs',         # 优化算法
    max_iter=1000,          # 最大迭代次数
    random_state=42,
    class_weight='balanced'  # 处理类别不平衡
)

# 6. 交叉验证评估模型
print("\n步骤6: 进行5折交叉验证")
# 使用分层K折交叉验证（保持类别分布）
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

# 评估多个指标
scoring_metrics = ['accuracy', 'precision', 'recall', 'f1', 'roc_auc']
cv_results = {metric: [] for metric in scoring_metrics}

for metric in scoring_metrics:
    scores = cross_val_score(
        model, X_train, y_train, 
        cv=cv, scoring=metric, n_jobs=-1
    )
    cv_results[metric] = scores
    print(f"  {metric.capitalize()}: 平均 {scores.mean():.3f} ± {scores.std():.3f}")

# 7. 训练最终模型
print("\n步骤7: 训练最终模型")
model.fit(X_train, y_train)

# 8. 在测试集上评估模型
print("\n步骤8: 在测试集上评估模型")
y_pred = model.predict(X_test)
y_prob = model.predict_proba(X_test)[:, 1]  # 预测生存概率

# 计算评估指标
accuracy = accuracy_score(y_test, y_pred)
roc_auc = roc_auc_score(y_test, y_prob)
precision = precision_score(y_test, y_pred)
recall = recall_score(y_test, y_pred)
f1 = f1_score(y_test, y_pred)

print(f"  测试集准确率: {accuracy:.4f}")
print(f"  ROC AUC: {roc_auc:.4f}")
print(f"  精确率: {precision:.4f}")
print(f"  召回率: {recall:.4f}")
print(f"  F1分数: {f1:.4f}")

# 9. 可视化分析
print("\n步骤9: 可视化分析")
plt.figure(figsize=(15, 10))

# 9.1 混淆矩阵
plt.subplot(2, 2, 1)
cm = confusion_matrix(y_test, y_pred)
plt.imshow(cm, interpolation='nearest', cmap=plt.cm.Blues)
plt.title('混淆矩阵')
plt.colorbar()
tick_marks = np.arange(2)
plt.xticks(tick_marks, ['死亡', '存活'])
plt.yticks(tick_marks, ['死亡', '存活'])
plt.xlabel('预测标签')
plt.ylabel('真实标签')

# 添加数值标签
thresh = cm.max() / 2.
for i in range(cm.shape[0]):
    for j in range(cm.shape[1]):
        plt.text(j, i, format(cm[i, j], 'd'),
                 horizontalalignment="center",
                 color="white" if cm[i, j] > thresh else "black")

# 9.2 ROC曲线
plt.subplot(2, 2, 2)
fpr, tpr, thresholds = roc_curve(y_test, y_prob)
plt.plot(fpr, tpr, label=f'Logistic Regression (AUC = {roc_auc:.2f})')
plt.plot([0, 1], [0, 1], 'k--')
plt.xlabel('假阳性率')
plt.ylabel('真阳性率')
plt.title('ROC曲线')
plt.legend(loc='lower right')

# 9.3 精确率-召回率曲线
plt.subplot(2, 2, 3)
precision_vals, recall_vals, _ = precision_recall_curve(y_test, y_prob)
ap = average_precision_score(y_test, y_prob)
plt.plot(recall_vals, precision_vals, label=f'Logistic Regression (AP = {ap:.2f})')
plt.xlabel('召回率')
plt.ylabel('精确率')
plt.title('精确率-召回率曲线')
plt.legend(loc='upper right')

# 9.4 特征重要性
plt.subplot(2, 2, 4)
coefficients = model.coef_[0]
feature_importance = pd.DataFrame({
    'Feature': selected_features,
    'Importance': np.abs(coefficients)
}).sort_values('Importance', ascending=True)  # 为水平条形图排序

# 创建水平条形图
y_pos = np.arange(len(feature_importance))
plt.barh(y_pos, feature_importance['Importance'], align='center', color='skyblue')
plt.yticks(y_pos, feature_importance['Feature'])
plt.xlabel('系数绝对值')
plt.title('特征重要性')

plt.tight_layout()
plt.savefig('horse_survival_analysis.png', dpi=300)
plt.show()

# 10. 模型解释
print("\n步骤10: 模型解释")
print("  逻辑回归系数 (正值表示增加生存概率):")
for feature, coef in zip(selected_features, coefficients):
    print(f"    {feature}: {coef:.4f}")

# 11. 模型部署建议
print("\n步骤11: 模型部署建议")
print("  基于分析结果，建议：")
print("  - 在临床环境中使用该模型作为辅助决策工具")
print("  - 重点关注关键特征：")
for feature in feature_importance.sort_values('Importance', ascending=False).head(5)['Feature']:
    print(f"      {feature}")
print("  - 定期使用新数据重新训练模型以保持性能")
print("  - 结合兽医专业知识解释模型预测结果")

print("\n分析完成！")